home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / csymm.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  9.7 KB  |  300 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE CSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB,
  5.      $                   BETA, C, LDC )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        SIDE, UPLO
  8.       INTEGER            M, N, LDA, LDB, LDC
  9.       COMPLEX            ALPHA, BETA
  10. *     .. Array Arguments ..
  11.       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  CSYMM  performs one of the matrix-matrix operations
  18. *
  19. *     C := alpha*A*B + beta*C,
  20. *
  21. *  or
  22. *
  23. *     C := alpha*B*A + beta*C,
  24. *
  25. *  where  alpha and beta are scalars, A is a symmetric matrix and  B and
  26. *  C are m by n matrices.
  27. *
  28. *  Parameters
  29. *  ==========
  30. *
  31. *  SIDE   - CHARACTER*1.
  32. *           On entry,  SIDE  specifies whether  the  symmetric matrix  A
  33. *           appears on the  left or right  in the  operation as follows:
  34. *
  35. *              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
  36. *
  37. *              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
  38. *
  39. *           Unchanged on exit.
  40. *
  41. *  UPLO   - CHARACTER*1.
  42. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  43. *           triangular  part  of  the  symmetric  matrix   A  is  to  be
  44. *           referenced as follows:
  45. *
  46. *              UPLO = 'U' or 'u'   Only the upper triangular part of the
  47. *                                  symmetric matrix is to be referenced.
  48. *
  49. *              UPLO = 'L' or 'l'   Only the lower triangular part of the
  50. *                                  symmetric matrix is to be referenced.
  51. *
  52. *           Unchanged on exit.
  53. *
  54. *  M      - INTEGER.
  55. *           On entry,  M  specifies the number of rows of the matrix  C.
  56. *           M  must be at least zero.
  57. *           Unchanged on exit.
  58. *
  59. *  N      - INTEGER.
  60. *           On entry, N specifies the number of columns of the matrix C.
  61. *           N  must be at least zero.
  62. *           Unchanged on exit.
  63. *
  64. *  ALPHA  - COMPLEX         .
  65. *           On entry, ALPHA specifies the scalar alpha.
  66. *           Unchanged on exit.
  67. *
  68. *  A      - COMPLEX          array of DIMENSION ( LDA, ka ), where ka is
  69. *           m  when  SIDE = 'L' or 'l'  and is n  otherwise.
  70. *           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
  71. *           the array  A  must contain the  symmetric matrix,  such that
  72. *           when  UPLO = 'U' or 'u', the leading m by m upper triangular
  73. *           part of the array  A  must contain the upper triangular part
  74. *           of the  symmetric matrix and the  strictly  lower triangular
  75. *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
  76. *           the leading  m by m  lower triangular part  of the  array  A
  77. *           must  contain  the  lower triangular part  of the  symmetric
  78. *           matrix and the  strictly upper triangular part of  A  is not
  79. *           referenced.
  80. *           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
  81. *           the array  A  must contain the  symmetric matrix,  such that
  82. *           when  UPLO = 'U' or 'u', the leading n by n upper triangular
  83. *           part of the array  A  must contain the upper triangular part
  84. *           of the  symmetric matrix and the  strictly  lower triangular
  85. *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
  86. *           the leading  n by n  lower triangular part  of the  array  A
  87. *           must  contain  the  lower triangular part  of the  symmetric
  88. *           matrix and the  strictly upper triangular part of  A  is not
  89. *           referenced.
  90. *           Unchanged on exit.
  91. *
  92. *  LDA    - INTEGER.
  93. *           On entry, LDA specifies the first dimension of A as declared
  94. *           in the  calling (sub) program. When  SIDE = 'L' or 'l'  then
  95. *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
  96. *           least max( 1, n ).
  97. *           Unchanged on exit.
  98. *
  99. *  B      - COMPLEX          array of DIMENSION ( LDB, n ).
  100. *           Before entry, the leading  m by n part of the array  B  must
  101. *           contain the matrix B.
  102. *           Unchanged on exit.
  103. *
  104. *  LDB    - INTEGER.
  105. *           On entry, LDB specifies the first dimension of B as declared
  106. *           in  the  calling  (sub)  program.   LDB  must  be  at  least
  107. *           max( 1, m ).
  108. *           Unchanged on exit.
  109. *
  110. *  BETA   - COMPLEX         .
  111. *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
  112. *           supplied as zero then C need not be set on input.
  113. *           Unchanged on exit.
  114. *
  115. *  C      - COMPLEX          array of DIMENSION ( LDC, n ).
  116. *           Before entry, the leading  m by n  part of the array  C must
  117. *           contain the matrix  C,  except when  beta  is zero, in which
  118. *           case C need not be set on entry.
  119. *           On exit, the array  C  is overwritten by the  m by n updated
  120. *           matrix.
  121. *
  122. *  LDC    - INTEGER.
  123. *           On entry, LDC specifies the first dimension of C as declared
  124. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  125. *           max( 1, m ).
  126. *           Unchanged on exit.
  127. *
  128. *
  129. *  Level 3 Blas routine.
  130. *
  131. *  -- Written on 8-February-1989.
  132. *     Jack Dongarra, Argonne National Laboratory.
  133. *     Iain Duff, AERE Harwell.
  134. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  135. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  136. *
  137. *
  138. *     .. External Functions ..
  139.       LOGICAL            LSAME
  140.       EXTERNAL           LSAME
  141. *     .. External Subroutines ..
  142.       EXTERNAL           XERBLA
  143. *     .. Intrinsic Functions ..
  144.       INTRINSIC          MAX
  145. *     .. Local Scalars ..
  146.       LOGICAL            UPPER
  147.       INTEGER            I, INFO, J, K, NROWA
  148.       COMPLEX            TEMP1, TEMP2
  149. *     .. Parameters ..
  150.       COMPLEX            ONE
  151.       PARAMETER        ( ONE  = ( 1.0E+0, 0.0E+0 ) )
  152.       COMPLEX            ZERO
  153.       PARAMETER        ( ZERO = ( 0.0E+0, 0.0E+0 ) )
  154. *     ..
  155. *     .. Executable Statements ..
  156. *
  157. *     Set NROWA as the number of rows of A.
  158. *
  159.       IF( LSAME( SIDE, 'L' ) )THEN
  160.          NROWA = M
  161.       ELSE
  162.          NROWA = N
  163.       END IF
  164.       UPPER = LSAME( UPLO, 'U' )
  165. *
  166. *     Test the input parameters.
  167. *
  168.       INFO = 0
  169.       IF(      ( .NOT.LSAME( SIDE, 'L' ) ).AND.
  170.      $         ( .NOT.LSAME( SIDE, 'R' ) )      )THEN
  171.          INFO = 1
  172.       ELSE IF( ( .NOT.UPPER              ).AND.
  173.      $         ( .NOT.LSAME( UPLO, 'L' ) )      )THEN
  174.          INFO = 2
  175.       ELSE IF( M  .LT.0               )THEN
  176.          INFO = 3
  177.       ELSE IF( N  .LT.0               )THEN
  178.          INFO = 4
  179.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  180.          INFO = 7
  181.       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
  182.          INFO = 9
  183.       ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
  184.          INFO = 12
  185.       END IF
  186.       IF( INFO.NE.0 )THEN
  187.          CALL XERBLA( 'CSYMM ', INFO )
  188.          RETURN
  189.       END IF
  190. *
  191. *     Quick return if possible.
  192. *
  193.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  194.      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  195.      $   RETURN
  196. *
  197. *     And when  alpha.eq.zero.
  198. *
  199.       IF( ALPHA.EQ.ZERO )THEN
  200.          IF( BETA.EQ.ZERO )THEN
  201.             DO 20, J = 1, N
  202.                DO 10, I = 1, M
  203.                   C( I, J ) = ZERO
  204.    10          CONTINUE
  205.    20       CONTINUE
  206.          ELSE
  207.             DO 40, J = 1, N
  208.                DO 30, I = 1, M
  209.                   C( I, J ) = BETA*C( I, J )
  210.    30          CONTINUE
  211.    40       CONTINUE
  212.          END IF
  213.          RETURN
  214.       END IF
  215. *
  216. *     Start the operations.
  217. *
  218.       IF( LSAME( SIDE, 'L' ) )THEN
  219. *
  220. *        Form  C := alpha*A*B + beta*C.
  221. *
  222.          IF( UPPER )THEN
  223.             DO 70, J = 1, N
  224.                DO 60, I = 1, M
  225.                   TEMP1 = ALPHA*B( I, J )
  226.                   TEMP2 = ZERO
  227.                   DO 50, K = 1, I - 1
  228.                      C( K, J ) = C( K, J ) + TEMP1    *A( K, I )
  229.                      TEMP2     = TEMP2     + B( K, J )*A( K, I )
  230.    50             CONTINUE
  231.                   IF( BETA.EQ.ZERO )THEN
  232.                      C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2
  233.                   ELSE
  234.                      C( I, J ) = BETA *C( I, J ) +
  235.      $                           TEMP1*A( I, I ) + ALPHA*TEMP2
  236.                   END IF
  237.    60          CONTINUE
  238.    70       CONTINUE
  239.          ELSE
  240.             DO 100, J = 1, N
  241.                DO 90, I = M, 1, -1
  242.                   TEMP1 = ALPHA*B( I, J )
  243.                   TEMP2 = ZERO
  244.                   DO 80, K = I + 1, M
  245.                      C( K, J ) = C( K, J ) + TEMP1    *A( K, I )
  246.                      TEMP2     = TEMP2     + B( K, J )*A( K, I )
  247.    80             CONTINUE
  248.                   IF( BETA.EQ.ZERO )THEN
  249.                      C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2
  250.                   ELSE
  251.                      C( I, J ) = BETA *C( I, J ) +
  252.      $                           TEMP1*A( I, I ) + ALPHA*TEMP2
  253.                   END IF
  254.    90          CONTINUE
  255.   100       CONTINUE
  256.          END IF
  257.       ELSE
  258. *
  259. *        Form  C := alpha*B*A + beta*C.
  260. *
  261.          DO 170, J = 1, N
  262.             TEMP1 = ALPHA*A( J, J )
  263.             IF( BETA.EQ.ZERO )THEN
  264.                DO 110, I = 1, M
  265.                   C( I, J ) = TEMP1*B( I, J )
  266.   110          CONTINUE
  267.             ELSE
  268.                DO 120, I = 1, M
  269.                   C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J )
  270.   120          CONTINUE
  271.             END IF
  272.             DO 140, K = 1, J - 1
  273.                IF( UPPER )THEN
  274.                   TEMP1 = ALPHA*A( K, J )
  275.                ELSE
  276.                   TEMP1 = ALPHA*A( J, K )
  277.                END IF
  278.                DO 130, I = 1, M
  279.                   C( I, J ) = C( I, J ) + TEMP1*B( I, K )
  280.   130          CONTINUE
  281.   140       CONTINUE
  282.             DO 160, K = J + 1, N
  283.                IF( UPPER )THEN
  284.                   TEMP1 = ALPHA*A( J, K )
  285.                ELSE
  286.                   TEMP1 = ALPHA*A( K, J )
  287.                END IF
  288.                DO 150, I = 1, M
  289.                   C( I, J ) = C( I, J ) + TEMP1*B( I, K )
  290.   150          CONTINUE
  291.   160       CONTINUE
  292.   170    CONTINUE
  293.       END IF
  294. *
  295.       RETURN
  296. *
  297. *     End of CSYMM .
  298. *
  299.       END
  300.